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The physics of an earthquake is a subject with many unknowns. It is true that 
we have a good understanding of the propagation of seismic waves through 
the Earth and that given a large set of seismographic records we are able 
to reconstruct a posteriori the history of the fault rupture (the origin of the 
waves). However, when we consider the physical processes which lead to the 
1-^ I initiation of a rupture with a subsequent slip and its growth through a fault 

r^ . system to give rise to an earthquake, then our knowledge is really limited. Not 

O I only the friction law and the rupture evolution rules are largely unknown, but 

the role of many other processes such as plasticity, fluid migration, chemical 
reactions, etc., and the couplings between them, remain unclear 0E]- 

On the other hand, one may wonder about the physics of many earth- 
quakes. How do the collective properties of the set defined by all earthquakes 
in a given region, or better, in the whole world, emerge from the physics of 
(y-v ' individual earthquakes? How does seismicity, which is the structure formed by 

■^ . all earthquakes, depend on its elementary constituents -the earthquakes? And 

^D [ which are these properties? Which kind of dynamical process does seismicity 

constitute? It may be that these collective properties are largely independent 
on the physics of the individual earthquakes, in the same way that many of 
t^ . the properties of a gas or a solid do not depend on the constitution of its 

elementary units -the atoms (for a broad range of temperatures it doesn't 
matter if we have atoms, with its complicated quantum structure, or mi- 
'^ ' croscopic marbles). It is natural then to consider that the physics of many 

C , earthquakes has to be studied with a different approach than the physics of 

one earthquake, and in this sense we can consider the use of statistical physics 
not only appropriate but necessary to understand the collective properties of 
earthquakes. 

^^ . Here, we provide a summary of recent work on the statistics of the tem- 

V-j ' poral properties of seismicity, considering the phenomenon as a whole and 

with the goal of looking for general laws. We show the fulfillment of a scaling 
law for recurrence-time distributions, which becomes universal for stationary 
seismicity and for aftershock sequences which are transformed into stationary 



196 Alvaro Corral 

processes by means of a nonlinear rescaling of time. The existence of a de- 
creasing power-law regime in the distributions has paradoxical consequences 
on the time evolution of the earthquake hazard and on the expected time of 
occurrence of an incoming event, as we will see. On the other hand, the scal- 
ing law for recurrence times is equivalent to the invariance of seismicity under 
renormalization-group-like transformations, for which the role of correlations 
between recurrence times and magnitudes is essential. Finally, we relate the 
recurrence-time densities studied here with the method previously introduced 
by Bak et al. |3]. 



1 The Gutenberg-Richter Law and the Omori Law 

Traditionally, the knowledge of seismicity has been limited to a few phe- 
nomenological laws, the most important being the Gutenberg-Richter (GR) 
law and the Omori law. The GR law determines that, for a certain region, 
the number of earthquakes in a long period of time decreases exponentially 
with the magnitude; to be concrete, N{Mc) ex 10^''*^'^, where N{Mc) is the 
number of earthquakes with magnitude M greater or equal than a threshold 
value Mc, and the 6— value is a constant usually close to one ^IHIEICI- 

If we introduce the seismic rate, r(i, Mc), defined as the number of earth- 
quakes with M > Mc per unit time in a time interval around t, then, the GR 
relation can be expressed in terms of the mean seismic rate, R{Mc), as 

RiMc) = (r(t, Mc)} = ^ / r{t, Mc)dt = ^^ = RolQ-'''^ (1) 



T Jo T 

where T is the total time under consideration and Rq is the (hypothetical) 
mean rate in the region for Mc = (its dependence, as well as that of other 
parameters, on the region selected for study is implicit and will not be indi- 
cated when it is superfluous). In fact, the GR law must be understood as a 
probabilistic law, and then we conclude that earthquake magnitude follows an 
exponential distribution, this is, Prob[A'f > Mc] = N{Mc)/M oc e-i"ioMfc^ 
with TV the total number of earthquakes, of any magnitude. Due to the prop- 
erties of the exponential distribution, the derivative of Prob[M > M^, which 
is the probability density (with a minus sign), is also an exponential. 

In terms of the seismic moment or of the dissipated energy, which are 
increasing exponential functions of the magnitude, the GR law transforms 
into a power-law distribution, the usual signature of scale invariance. This 
means that earthquakes have no characteristic size of occurrence, if we take 
the seismic moment or the energy as more appropriate measures of earthquake 
size than the magnitude 0j. 

The Omori law (in its modified form) states that after a strong earthquake, 
which is called mainshock, the seismic rate for events with M > Ale in a 
certain region around the mainshock increases abruptly and then decays in 
time essentially as a power law; more precisely. 
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K^,M.) = /^^ (2) 

^ ' ^ (1 + t/c)P ^ ' 

where t is the time measured from the mainshock, ro(Afc) is the maximum rate 
for M > Mci which coincides with the rate immediately after the mainshock, 
i.e., r(i = 0, Mc) = ro(Mc), c is a short-time constant (of the order of hours or 
a few days), which describes the deviation from a pure power law right after 
the mainshock, and the exponent p is usually close to 1. In fact, c depends to 
a certain degree on Mc and, together with ro(Mc) and p, depends also on the 
mainshock magnitude [315117]. 

Nowadays it has been confirmed that the Omori law does not only apply to 
strong earthquakes, but to any earthquake, with a productivity factor (rg) for 
small earthquakes which is orders of magnitude smaller than for large events. 
In this way, the classification of earthquakes in mainshocks and aftershocks 
turns out to be only relative, as we will have a cascade process in which after- 
shocks become also mainshocks of secondary sequences and so on. When an 
aftershock happens to have a magnitude larger than the mainshock a change 
of roles occur: the mainshock is considered a foreshock and the aftershock 
becomes the mainshock. Also, the triggering of strong aftershocks may cause 
that the overall seismic rate departs significantly from the Omori law, as it 
happens in earthquake swarms. 

In any case, the Omori law illustrates clearly the temporal clustering of 
earthquakes, for which events (aftershocks) tend to gather close (in time) 
to a strong event (the mainshock), becoming more dilute as time from the 
mainshocks grows. In addition, the fact that the seismic rate decays essentially 
as a power law means that the relaxation process has no characteristic time, in 
opposition to the usual situation in physics (think for instance in radioactive 
decay). Finally, the Omori law has a probabilistic interpretation, as an all- 
return-time distribution, measuring the probability that earthquakes occur at 
a time t after a mainshock. 



2 Recurrence-Time Distributions and Scaling Laws 

One can go beyond the GR law and the Omori law and wonder about the tem- 
poral properties of individual earthquakes (from a statistical point of view), 
in particular about the time interval between consecutive earthquakes. In this 
case, it is necessary to assume that earthquakes are point events in time, or at 
least that their temporal properties are well described by their initiation time. 
In contrast to the previous approaches, this perspective has been much less 
studied and no general law has been proposed; rather, the situation is con- 
fusing in the literature, where claims range from nearly-periodic behavior for 
large earthquakes to totally random occurrence of mainshocks (see the cita- 
tions at Refs. |10[lll| '). Furthermore, it can be argued that the times between 
consecutive earthquakes depend strongly on the selection of the coordinates 
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of the region under study and the range of magnitudes selected (which change 
the sequence of events) and therefore one is dealing with an ill-defined vari- 
able. We will see that the existence of universal properties for these times 
invalidates this objection. 

Following the point of view of Bak et al. , we have addressed this problem 
by considering seismicity as a phenomenon on its own. In this way, we will 
not separate events into different kinds (foreshocks, mainshocks, aftershocks, 
or microearthquakes, etc.), nor divide the crust into provinces with different 
tectonic properties, but will place all events and regions on the same footing; 
in other words, we wonder about the very nature of seismicity as a whole, 
from a complex-system perspective, in opposition to a reductionist approach 

This exposition will concentrate on the temporal properties of seismicity, 
and their dependence with space and magnitude, but equally important are 
the spatial properties. It turns out that all the aspects of seismicity are closely 
related to each other and one cannot study them separately. Although all the 
events are important, as they are the elementary constituents of seismicity, we 
will need to consider windows of observation in space, time, and magnitude; of 
course, this is due to the incompleteness of seismic records but also to the fact 
that the variation of the quantities we measure with the range of magnitudes 
selected or with the size of the spatial region under study will allow us to 
establish self-similar properties for seismicity. 

Let us select an arbitrary region of the Earth, a temporal period, and a 
minimum magnitude Mc, in such a way that only events in this space-time- 
magnitude window are taken into account. We can consider the resulting 
events as a point process in time, disregarding the magnitude and the spatial 
degrees of freedom (this is not arbitrary, as Mc and the size of the region 
will be systematically varied later on), in this way we can order the events in 
time, from i = 1 to N{Mc) and characterize each one only by its occurrence 
time, ti. From here we can define the recurrence time r (also called waiting 
time, interevent time, interoccurrence time, etc.) as the time interval between 
consecutive events, i.e., Ti = ti — ti^i. The mean recurrence time, {t{Mc)), is 
obviously given by the inverse of the rate, i?~^(Mc); however, as the recur- 
rence time is broadly distributed, the mean alone is a poor characterization 
of the process and it is inevitable to work with the probability distribution of 
recurrence times. So, we compute the recurrence-time probability density as 

„ , , Probfr < recurrence time < r + dr] , , 

D{t; Mc) - ^ — = ^, (3) 

where dr has to be small enough to allow D to represent a continuous function 
but large enough to contain enough data to be statistically significant (note 
that the spatial dependence of D is not indicated explicitly) . 
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2.1 Scaling Laws for Recurrence- Time Distributions 

We can illustrate this procedure with the waveform cross-correlation catalog 
of Southern California obtained by Shearer et al.^ 13. for the years 1984- 
2002, containing 26700 events with M > 2.5 (84209 events with M > 2). 
The recurrence-time probability densities for several values of Mc are shown 
in Fig. ^left). First, one can see that r ranges from seconds to more than 
100 days (in fact, we have restricted our analysis to recurrence-times greater 
than one minute; shorter times do not follow the same trend than the rest, 
probably due to the incompleteness of the records in that time scale). Also, 
the different distributions look very similar in shape, although the ranges are 
different (obviously, the larger Mc, the smaller the number of events N(Mc), 
and the larger the mean time between them). 




T (mill) 




Fig. 1. (left) Probability densities of recurrence times in Southern California (SC) 
for the period 1984-2002, for several Mc values, (right) The same probability densi- 
ties rescaled by their rate. The data collapse illustrates the fulfillment of a scaling 
law. The continuous line is a gamma fit. 



Figure fright) shows the same distributions but rescaled by the mean 
rate, as a function of the rescaled recurrence time, i.e., D{t; Mc)/R{Mc) versus 
R{Mc)t. In this case all the distributions collapse onto a single curve / and 
we can establish the fulfillment of a scaling law |14| , 



D{t; Mc) = R[Mc)f{R{Mc)T) 



(4) 



where / is the scaling function, and corresponds to the recurrence-time density 
in the hypothetical case R{Mc) — 1. Note that we could have arrived to a 
similar equation by scaling arguments, but there would be no reason for the 
function / to be independent on Mc- Only imposing the self-similarity of the 
process in time-magnitude can lead to the fact that / does not depend on Mc 
and therefore to the fact that / is a scaling function. As R{Mc) verifies the 
GR law, the scaling law can be written 



^Available at |http : //-mra . data . scec . org/f tp/catalogs/SHLK/ 1 
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D{t; M,) = 10-''^^=/(10-''^^^r). (5) 

The GR law can be calculated from the scaling law; just calculate the mean 
recurrence time, {T{h'Q) = J^^ tD{t; AQdr = 10*^= /(,°° z/(z)dz oc 10^^^ 
and as {t{Mc)) is the inverse of the mean rate, then, R{Mc) oc 10"^^'=. But 
the scaling law does not only include the GR law, it goes one step further, 
as it implies that the GR law is fulfilled at any time, if times are properly 
selected; indeed, events separated by recurrence times r' for M > M^ and r 
for M >Mc occur at a GR ratio, 10-^(*^o-Mc)^ jf ^nd only if the ratio of the 
recurrence times is given by 10^(*^c^^^=). Notice that the only requirement for 
the GR law to be fulfilled (for a long period of time) is that D{t; Mc) has a 
mean that verifies the GR law, i.e., {t{Mc)) = R~^{Mc) = RoW'^''^; therefore, 
the fulfillment of the GR law at any time is a new feature of seismicity. 

To make it more concrete, we can count the number of events in Southern 
California with M > 3 coming after a recurrence time r = 100 hours and 
compare with the number of events with M > A after the same recurrence 
time; then the ratio of these numbers has nothing to do with the GR relation. 
However, if for M > 4 we select events with r = 1000 hours (the 6— value in 
the GR law is very close to 1 in Southern California) then, the number of these 
events is about 1/10 of the number of events with M > 3 and r — 100 hours, 
the same proportion as when we consider all events (no matter the value of 
r). This could be somehow analogous to the well-known law of correspond- 
ing states in condensed-matter physics: two pairs of consecutive earthquakes 
in different magnitude windows would be in "corresponding states" if their 
rescaled recurrence times are the same. 

2.2 Relation with the Omori Law 

In general, as seismicity is not stationary, the scaling function / will change 
with the spatio-temporal window of observation. In the case of Omori after- 
shock sequences, the scaling function, and therefore the distribution of recur- 
rence times, is related to the Omori law, as we now see. Let us assume, just for 
simplicity, that the aftershock sequence can be modeled as a nonhomogeneous 
Poisson process (also called nonstationary Poisson process, this is a Poisson 
process but with a time-variable rate, in such a way that at any instant the 
probability of occurrence, per unit time, is not constant but is independent on 
the occurrence of other events) ; in this case the rate of occurrence will be given 
by the Omori law, Eq. (j^J. Then, the recurrence-time density is a temporal 
mixture of Poisson processes, which have a density D{T\r{Mc)) = re^*"^, so, 

1 /""" 1 r° 

D{t-M^) = -J rD{T\r)p{r;Mc)dr= - r'^e'"'^ p{r; Mc)dr, (6) 

where p(r; Mc) is the density of rates, /i is a normalization factor that turns 
out to be the mean value oi r, fj, = {r{Mc)) = J rp{Mc)dr and ro(Mc) and 
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TmiMc) the maximuni and minimuni rate, respectively, assuming vq 3> r,„; 
the factor r appears because the probabiUty of a given Z)(r|r) to contribute 
to D{t;Mc) is proportional to r. 

The density of rates can be obtained by the projection of r{t; Mc) onto the 
r axis, turning out to be. 



p(r; M,) ex 



dr 



C 
p{r; Mc) = ^T^T^ for r,n < r < Vq (7) 



with C just a constant (depending on Mc) that can be obtained from normal- 
ization. Substituting, we get 



(~i fro 



C[r{2 - l/p, r„,T) - r{2 - 1/p, ror)] 



//t2-1/p 



(8) 
with r(a,z) = J z" ^e ^dz the incomplete gamma function (note that 
r{l,z) = e~^). It is clear that for intermediate recurrence times, l/rp <^ r ^ 
l/j'm, we get a power law of exponent 2 — l/p for the recurrence time density, 

with r{a) the usual (complete) gamma function. This power-law behavior 
has been derived before by Senshu and by Utsu for nonhomogeneous Poisson 
processes 0, but our procedure can be easily extended beyond this case, just 
defining a different D(T\r), for which the value of the recurrence-time exponent 
2 — l/p is still valid. Notice that the value of this exponent is close to one if 
the p— value is close to one, but both exponents are only equal if p = 1, in any 
other case we have 2 — l/p < p, which means that in general D{t; Mc) decays 
more slowly than r(t; Mc). If we consider large recurrence times, rmT 3> 1, we 
can use the asymptotic expansion r{a, z) — > z"~^e~^ + ■ ■ ■ for z -^ oo jl5|. 
to get 

D{r;Mc)^^^^^^^^- , (10) 

which in the limit we are working is essentially an exponential decay. 

Although the equations derived here for a nonhomogeneous Poisson pro- 
cess with Omori rate reproduce well the recurrence-time distribution of after- 
shock sequences, D{t; Mc) [I2j, the choice of an exponential form for Z?(r|r) is 
not justified, as we will see in the next sections. Nevertheless, for the moment 
we are only interested in the form of D{t; Mc). 

2.3 Gamma Fit of the Scaling Function 

The fact that the density D{t; Mc) for a nonhomogeneous Poisson-Omori 
sequence is a power law for intermediate times and follows Eq. Hl()|l for long 
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times suggests that a simpler parameterization of the distribution can be 
obtained by the combination of both behaviors; in the case of the scaling 
function /, which must follow the same distribution as D (but with mean 
equal to one), we can write 

p-e/a 
f{0) « 



6'2-i/p(l + 6i)i/p-i' 

However, as both power laws of 9 are very similar and in the long time limit it 
is the exponential alone what is really important, we can simplify even further 
and use the gamma distribution to model /; so, 

where 6 plays the role of a dimensionless recurrence time, 9 = Rt, a is a 
dimensionless scale parameter, and C is a correction to normalization due to 
the fact that the gamma distribution may not be valid for very short times; 
this will allow the shape parameter 7 not to be restricted to the case 7 > 0, 
the usual condition for the gamma distribution (nevertheless, if 7 < the 
factor 1^(7) is inappropriate for normalization). As / is introduced in such a 
way that the mean of 9 is {9) = 1, the parameters are not independent; for 
instance, for C = 1, (9) = ja and in consequence a = I/7. So, essentially, we 
only have one parameter to fit, 7, to characterize the process. In the case of 
Omori sequences, 1 — 7 = 2 — 1/p and a = R/rm, ^ 7 ~ Tm/R, but we will 
see that the gamma distribution has a wider applicability than just Omori 
sequences. 

A fit of the gamma distribution to the rescaled distribution for Southern- 
California, shown in Fig. fright), yields the parameter values 7 ~ 0.22 and 
a ~ 3; this yields a power-law exponent for small and intermediate times 
1 — 7 ~ 0.78 and allows to calculate a p— value p — {1 + 7)^^ — 0.82, which 
can be interpreted as an average for Southern California, and a minimum rate 
r„i — R/3. Of course, with our resolution we only can establish 1— 7 ~ p ~ 0.8. 

2.4 Universal Scaling Law for Stationary Seismicity 

We have mentioned the nonstationary character of seismicity and that in 
consequence the scaling function / depends on the window of observation. 
A more robust, universal law can be established if we restrict our study to 
stationary seismicity. By stationary seismicity we mean in fact homogeneity 
in time, which implies that the statistical properties of the process do not 
depend on the time window of observation, in particular, the mean rate must 
be practically constant in time. 

It is obvious that an aftershock sequence following the Omori law (with 
p > 0) is not stationary, but observational evidence shows that in other cases 
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seismicity can be well described by a stationary process, for example world- 
wide seismicity for the last 30 years (for which there are reasonably good 
data) or regional seismicity in between large aftershock sequences. It should 
be clear that considering stationary seismicity has nothing to do with declus- 
tering (the removal of aftershocks from data). We simply consider periods of 
time for which no aftershock sequence dominates in the spatial region selected 
for study, but many smaller sequences may be hidden in the data, intertwined 
in such a way to give rise to an overall stationary seismic rate. 

The total number of earthquakes in Southern-California (from Shearer 
et al.'s catalog) as a function of time since 1984 is displayed in Fig. |21 
Clearly, the behavior of the number of earthquakes in time is nonlinear, 
with episodic abrupt increments which correspond to large aftershock se- 
quences, following the trend prescribed by the Omori law, N{Mc,t) — 
N{Mc,0) + L ro(Mc)/(l +t' /c)Pdt'. However, there exist some periods which 
follow a linear increase of N{Mc,t) versus t; in particular, we have cho- 
sen for analysis the intervals (in years, with decimal notation) 1984-1986.5, 
1990.3-1992.1, 1994.6-1995.6, 1996.1-1996.5, 1997-1997.6, 1997.75-1998.15, 
1998.25-1999.35, 2000.55-2000.8, 2000.9-2001.25, 2001.6-2002, and 2002.5- 
2003. These intervals comprise a total time span of 9.25 years and contain 
6072 events for M > 2.5, corresponding to a mean rate i?(2.5) — 1.7 earth- 
quakes/day. Note from the figure that not only the rate of occurrence is nearly 
constant for each interval, but different intervals have similar values of the 
rate. 
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Fig. 2. Accumulated number of earthquakes in Southern California as a function 
of time. Some stationary or nearly stationary periods mentioned in the text are 
specially marked, see subsection 15.81 
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We will study all these stationary periods together, in order to improve the 
statistics. The probability densities of the recurrence times are calculated from 
all the periods and the corresponding rescaled distributions appear in Fig. 
ISJleft). The good quality of the data collapse indicates the validity of a scaling 
law of the type of Eq. |@J, although the scaling function / is clearly different 
than the one for the whole time period 1984-2002, in particular, the power-law 
is much flatter, which is an indication that the clustering degree is smaller in 
this case, in comparison, but still exists. Nevertheless, it is remarkable that 
this kind of clustering is different than the clustering of aftershock sequences, 
as in this case we are dealing with a stationary process. The figure shows also 
a plot of the scaling function / parameterized with a gamma distribution with 
7 = 0.7 and a — 1.38, which indeed implies a power-law exponent 1 — 7 = 0.3. 
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Fig. 3. Rescaled recurrence-time probability densities for the stationary periods 
explained in the text for Southern California (left) and for worldwide seismicity 
(right). The solid line is the same function in both cases, showing the universal 
character of the scaling law fulfilled. 



We now present the results for recurrence times in worldwide scale, using 
the NEIC-PDE worldwide catalog (National Earthquake Information Center, 
Preliminary Determination of Epicenters ^) which covers the period 1973-2002 
and yields 46055 events with M > 5. In this case the total number of earth- 
quakes grows linearly in time, which confirms the stationarity of worldwide 
seismicity. The corresponding rescaled recurrence-time probability densities 
are shown in Fig.|2Iright), together with the scaling function used in the pre- 
vious case (i.e., Southern-California stationary seismicity). The collapse of the 
data onto a single curve is again an indication of the validity of a scaling law, 
and the fact that this curve is well fit by the same scaling function than in the 
Southern-California stationary case is a sign of universality. We use the term 
universality with the usual meaning in statistical physics, in which it refers to 
very different systems (gases or magnetic solids, or in our case seismic occur- 



^ Available at |http : //wwwneic . or . usgs . gov/neis/epic/epic_global . html| 
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rence in quite diverse tectonic environments) sharing the sanre quantitative 
properties. 

In fact, the universality of the scaling law for recurrence-time distributions 
in the stationary case has been tested for several other regions, namely, Japan, 
Spain, New Zealand, New Madrid (USA), and Great Britain, with magnitude 
values ranging from Af > 1.5 to M > 7.5 (which is equivalent to a factor 10^ in 
the minimum dissipated energy), and for spatial areas as small as 0.16° ~ 20 
km [Hini 

2.5 Universal Scaling Law for Omori Sequences 

We now return to nonstationary seismicity to show how the universal scaling 
law for recurrence times applies there. For this purpose, let us consider the 
Landers earthquake, with magnitude M — 7.3, the largest event in Southern 
California in the last decades, taking place in 1992, June 28, at 34.12°A^, 
116.26°M^. After the earthquake, seismicity in Southern California followed 
the usual behavior when large shallow events happen: a sudden enormous 
increase in the number of earthquakes and a consequent slow decay in time, 
in good agreement with the Omori law. 

The previous universal results for stationary seismicity can be generalized 
in the nonstationary case by replacing the mean seismic rate R{Mc) by the 
"instantaneous" seismic rate r(i, M^) as the scaling factor in Eq. Q. Then, in 
order to obtain the rescaled, dimensionless recurrence time 0, we will rescale 
each recurrence time as 

e,=r{U-M,)n. (12) 

This means that it is the instantaneous rate of occurrence which sets the time 
scale. 

First, we examine the complete seismicity with M > Mc for a square region 
(in a space in which longitude and latitude are considered as rectangular 
coordinates), the region containing the Landers event (but not centered on 
it); that is, as in previous sections, we will not separate aftershocks from the 
rest of events. As expected, after a few days from the mainshock, the seismic 
rate decreases as a pure power law, which is equivalent to take i 3> c in Eq. 
(O, so, 

/C\P 

r(i;Afe)-ro(-j , 

which lasts until the rate reaches the background seismic level. This form for 
r{t;Mc) is fit to the measured seismic rate, see Fig.0Jright). One advantage 
of analyzing the pure power-law regime only, rather than the whole sequence 
using the modified Omori law, Eq. (O, is that it is believed that the deviations 
from power-law behavior for short times are due to the incompleteness of the 
catalogs after strong events; therefore, in this way we avoid the problem of 
incompleteness. 

Next, using the results of the fit rather than the direct measurement of 
r{t; Mc) we calculate 
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In fact, this rescaling could be replaced by 9i — r{ti- 
Ti/2)Ti, with no noticeable difference in the results, 
slowly at the scale of the recurrence time. 



(13) 

i)Ti or by Oi = r(ti_i + 
as the rate varies very 
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Fig. 4. (left) Seismic rate as a function of the time elapsed since the Landers 
earthquake for regions of different size L including the event, using the SCSN catalog. 
Only events with M > 2 are considered. The straight lines correspond to power laws, 
with exponents given by the p— value, (right) Recurrence-time probability densities 
for the power-law regime of the decay of the rate after the Landers event, rescaled 
at each time by the rate. The solid line represents the universal scaling function in 
terms of a gamma distribution. 



The probability densities of the rescaled recurrence times 9i obtained in 
this way are displayed in Fig. fright), showing a slow power-law decay fol- 
lowed by a faster decay, in surprising agreement, not only qualitative but also 
quantitative, with the results for stationary seismicity, in such a way that 
the universal scaling function for the stationary case is still valid |14l nT)] . 
Therefore, as in that case, the power-law regime in the density is a sign of 
clustering, but as the primary clustering structure of the sequence has been 
removed by the rescaling with Oniori rate, this implies the existence of a sec- 
ondary clustering structure inside the main sequence, due to the fact that any 
large aftershock may generate its own aftershocks [181 [TO] . What is remark- 
able is that this structure seems to be identical to the one corresponding to 
stationary seismicity. 

An important consequence of this is that the time behavior of seismicity 
depends on just one variable: the seismic rate. Another implication is the fact 
that aftershock sequences cannot be described as a nonhomogeneous Poisson 
process, as in that case one should obtain an exponential distribution for f{9). 
The use of the nonhomegeneous Poisson process previously in this work must 
be understood only as a first approximation to justify the use of the gamma fit 
for the distribution of recurrence times in an Omori sequence. Nevertheless, 
we will see in the next sections that the generalization of the nonhomogeneous 
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Poisson process taking into account the results explained here leads to similar 
conclusions for the time distribution in the sequence. 

The rescaling of the recurrence times with the seismic rate r{t, Ale) can 
be applied also to the occurrence times ti, in order to transform the Omori 
sequence (or in general any sequence with a time- variable rate) into a station- 
ary sequence. For this purpose we define the accumulated rescaled recurrence 
time 0, defined as 6*^ = 9i + 62 + ■ ■ ■ + Oi, which plays the role of a station- 
ary occurrence time, in the same way that 9i plays the role of a stationary 
recurrence time, in general. This allows the complete comparison between 
stationary seismicity and aftershock sequences. 



3 The Paradox of the Decreasing Hazard Rate and the 
Increasing Time Until the Next Earthquake 

Other functions, in addition to the probability density, are suitable for describ- 
ing the general properties of recurrence times. Although from a mathematical 
point of view the functions we are going to introduce are fully equivalent to 
the probability density, they show much clearly some interesting temporal 
features of seismicity. 

3.1 Decreasing of the Hazard Rate 

Let us consider first the hazard rate, A(t;Mc), defined for a certain region 
and for M > Mc as the probability per unit time of an immediate earthquake 
given that there has been a period r without activity [20] , 

w ^^x P£^b[r<T'<r + dr|T'>r] _D{t;M,) 



dT S{t; Mc) ' 

where r' is a generic label for the recurrence time, while r refers to a par- 
ticular value of the same quantity, the symbol | denotes conditional prob- 
ability, and S{t; Mc) is the survivor function, 5(t; M^) = Prob[T' > r] = 
J D{T'\Mc)dT' . Introducing the scaling law Q for D in the definitions it 
is immediate to obtain that both S{t; M^) and A(r; M^) verify also scaling 
relations, S{t;Mc) = g{RT) and X{t;Mc) — Rhi^Rr). If we make use of the 
gamma parameterization, Eq. Hll|l. with C ~ 1, we get for the scaling function 
/i, 

h{e) = - ^^ 



rii,o/a)- 

For short recurrence times this function diverges as a power law 

^ ' ~ r{-f)a'>e^''r 
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(in fact, in this limit the hazard rate becomes undistinguishable from the 
probabihty density). On the other hand, h{9) tends as a power law to the value 
1/a as -^ cxd; indeed, making use of the expansion r{'-f,z) -^ z^~^e~^[\ — 
(1 - 7)/z + (1 - 7)(2 - 7)/z2 + . . .] [TSl, we get 



h{e) = - 

a 



a(l-7) 
1+ , ^^ 



The overall behavior for 7 < 1 is that h{6) decreases monotonically as 9 
increases; so, contrary to common belief and certainly counterintuitively, these 
calculations allow us to predict that the hazard does not increase with the 
elapsed time since the last earthquake, but just the opposite, it decreases up 
to an asymptotic value that corresponds to a Poisson process of rate R/ a. 
This means that although the hazard rate decreases, it never reaches the zero 
value, and sooner or later a new earthquake will strike. 

If we compare the hazard rate for 7 < 1 with that of a Poisson process 
with the same mean (given by 7 = a = 1, and rate R), we see that for short 
recurrence times the hazard rate is well above the Poisson value, implying that 
at any instant the probability of having an earthquake is higher than in the 
Poisson case. In contrast, for long times the probability is below the Poisson 
value, by a factor 1/a. This is precisely the most direct characterization of 
clustering in time, for which we can say that events tend to attract each 
other, being closer in short time scales and more separated in long time scales 
(in comparison with the Poisson process). 

In conclusion, we predict that seismicity is clustered independently on the 
scale of observation; in the case of stationary seismicity this clustering is much 
less trivial than the clustering due to the increasing of the rate in aftershock 
sequences. Taking advantage of the self-similarity implied by the scaling law, 
we could extrapolate the clustering behavior to the largest events worldwide 
(M > 7.5) over relatively small spatial scales (hundreds of kilometers) and we 
would obtain a behavior akin to the long-term clustering observed by other 
means |2J 

3.2 Increasing of the Residual Time Until the Next Earthquake 

Let us introduce now the expected residual recurrence time, e{To;Mc), which 
provides the expected time till the next earthquake, given that a period tq 
without earthquakes (in the spatial area and range of magnitudes considered) 
has elapsed [23], 



e(To; M,) = (r - To I T > To) - — — - / (r - to)D{t; M,)dT. 

S{to;Mc) Jro 

where | denotes that the mean is calculated only when the condition r > tq is 
fulfilled. Again, the scaling law for D implies a scaling form for this function, 
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which is e(To) — e{RT^)) / R, and introducing the gamma parameterization we 
get for the scahng function 



r(7+l,g/a) 



7+ - 



7 Q-e/a 

a) r{^,e/a) 



making use of the relation 1^(7 + 1, z) = 7/^(7, z) + z'^'e ^. For short times we 
obtain 

remember that the unconditional mean is {9) = 7a = 1, precisely the value 
obtained for 9 = 0. For long times e{9) reaches, again as a power law, an 
asymptotic value equal to a, i.e., 



e{9) = a 



^ «(l-7) 



The global behavior of e{9) is monotonically increasing as a function of 9 if 
7 < 1. Therefore, the residual time until the next earthquake should grow with 
the elapsed time since the last one. Notice the counterintuitive behavior that 
this represents: if we decompose the recurrence time t as t = tq + tj, with Tf 
the residual time to the next event, the increase of tq implies the increase of 
the mean value of tj, but the mean value of r is kept fixed. In fact, this is fully 
equivalent to the previously reported decreasing-hazard phenomenon and just 
a more dramatic version of the classical waiting-time paradox |22[ I23L 1^ . 

This result seems indeed paradoxical for any time process, as we natu- 
rally expect that the residual recurrence (or waiting) time decreases as time 
increases; think for instance that you are waiting for the subway: you are con- 
fident that the next train is approaching; or when you celebrate your birthday, 
your expected residual lifetime decreases (at any time, in fact). Of course, for 
an expert statistician the case of earthquakes is not paradoxical, but only 
counterintuitive, and he or she can provide the counterexamples of newborns 
(mainly in underdeveloped countries) or of private companies, which become 
healthier or more solid as time passes and therefore their expected residual 
lifetime increases with time. These counterintuitive behaviors can be referred 
to as a phenomenon of negative aging. 

Nevertheless, for the concrete case of earthquakes the increasing of the 
expected residual recurrence time is still paradoxical, since one naively ex- 
pects that the longer the time one has been waiting for an earthquake, the 
closer it will be, due to the fact that as time passes stress increases on the 
faults and the next earthquake becomes more likely. Nevertheless, note that 
our approach does not deal with individual faults but with two-dimensional, 
extended regions, and in this case the evolution of the stress is not so clear. It 
is worth mentioning that, as far as the author knows, no conclusive study of 
this kind has been performed for observational data in individual faults, the 
difficulty on associating earthquakes to faults is one the major problems here. 
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3.3 Direct Empirical Evidence 

Our predictions for earthquake recurrence times follow the line initiated by 
other authors. Davis et al. |5S], pointed out that when a lognormal distribu- 
tion is a priori assumed for the recurrence times, the expected residual time 
increases with the elapsed time. However, the increase there was associated 
to the update of the distribution parameters as the time since the last earth- 
quake (which was taken into account in the estimation) increased, and not 
to an intrinsic property of the distribution. Sornette and Knopoff |26| showed 
that the increase (or decrease) depends completely on the election of the dis- 
tribution, and studied the properties of a number of them. We now will see 
that the observational data provide direct and clear evidence in favor of the 
picture of an incoming earthquake which is moving away in time. 
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Fig. 5. Rescaled hazard rate (left) and rescaled expected residual recurrence time 
(right) as a function of time for Southern California, 1988-1991 (nearly stationary 
period) and for worldwide seismicity, 1973-2002. The observational data agrees with 
the scaling functions derived from the gamma distribution. In the right plot the 
parameter a is not free, but a — I/7 to enforce e(0) = 1. 



Indeed, in order to rule out the possibility that these paradoxical predic- 
tions are an artifact introduced by the gamma parameterization, we must con- 
trast them with real seismicity; in fact, both the hazard rate and the expected 
residual recurrence time can be directly measured from the catalogs, with no 
assumption about their functional form. Their definitions provide a simple 
way to estimate these functions, and in this way we have applied these defini- 
tions to the recurrence-time data ^^I ■ From the results displayed in Fig. it 
is apparent that in all cases the hazard rate decreases with time whereas the 
expected residual recurrence time increases, as we have predicted. Although 
both quantities are well approximated by the proposed universal scaling func- 
tions, we emphasize that their behavior does not depend on any modeling of 
the process and in particular is independent on the gamma parameterization. 
Moreover, the fact that e(To) is far from being constant at large times means 
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that the time evolution is not properly described by a Poisson process, even 
in the long-time Hmit. 

We conclude stating that the contents of this section can be summarized in 
this simple sentence: the longer since the last earthquake, the lowest the hazard 
for a new one, which is fully equivalent to this one (although less shocking): 
the longer since the last earthquake, the longer the expected time till the next. 
Moreover, this happens in a self-similar way, thanks to the scaling laws which 
are fulfilled. 



4 Scaling Law Fulfillment as Invariance Under a 
Renormalizat ion- Group Transformation 

It is interesting to realize that the scaling law for the recurrence-time distribu- 
tion, Eq. 0), implies the invariance of the distribution under a renormalization- 
group transformation. Let us investigate deeper the meaning of the scaling 
analysis we have performed and its relation with the renormalization group. 




loon 1LILIII2 I'ICJO.I I'J'JIin I'J'JIIS 1991 

time (years A.D. J 



1990 1990.2 1990.1 1990.6 1990.8 1991 

time (years A.D.) 




I'i'ii 1 

time (\eals A D ) 



Fig. 6. Magnitude versus time of occurrence of worldwide earthquakes for sev- 
eral magnitude-time windows. The rising of the magnitude threshold from 5 to 6 
illustrates the thinning or decimation process characteristic of the first step of a 
renormalization-group transformation. The second step is given by the extension 
(rescaling) of the time axis from one year (1990) to 10 years (1990-2000). Notice the 
similarity between the first plot and the last one, which is due to the invariance of 
seismicity under this transformation. 



Figure |^ displays the magnitude M versus the occurrence time t of all 
worldwide earthquakes with M > Mc for different periods of time and Mc 
values. The top of the figure is for earthquakes with Af > 5 for the year 
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1990. If we rise the threshold up to M^ = 6 we get the results shown in 
Fig. Ef medium) . Obviously, as there are less earthquakes in this case, the 
distribution of recurrence times (time interval between consecutive "spikes" 
in the plot) becomes broader with respect to the previous case, as we know. 
The rising of the threshold can be viewed as a mathematical transformation 
of the seismicity point process, which is referred to as thinning in the context 
of stochastic processes |27| and is also equivalent to the common decima- 
tion performed for spin systems in renormalization-group transformations 
| 28l 1291 150] ; the term decimation is indeed appropriate as only one tenth of 
the events survive this transformation, due to the fulfillment of the GR law 
with 6=1. Figure inibottom) shows the same as Fig. |HI medium) but for ten 
years, 1990-1999, and represents a scale transformation of seismicity (also as 
in the renormalization group), contracting the time axis by a factor 10 to 
compensate for the previous decimation. The similarity between Fig. El^top) 
andl^bottom) is apparent, and is confirmed when the probability densities of 
the corresponding recurrence times are calculated and rescaled following Eq. 
(@J, see again Fig.lSfright). 

4.1 Simple Model to Renormalize 

A simple model may illustrate these ideas |31]. Let us assume that seismicity 
could be described as a time process for which each recurrence time r^ (which 
separates event i — 1 and i) only depends on Mi^i, the magnitude of the last 
event that has occurred before event i. Any other dependences are ignored, 
and in particular the values of the magnitudes are generated independently 
from the rest of the process. It is possible to shown that for this process the 
recurrence-time density for events above M^, D{t\ M^), can be related to the 
recurrence-time density for events above Mc conditioned to Mpre > M'^ or to 
Mpre < M^, which we denote D{T\Mpre > M^; AQ and D{T\Mpre < M'^; M^), 
respectively, where Mpre refers to the magnitude of the event immediately 
previous to the recurrence time, and it is assumed that M'^ > M^- The relation 
turns out to be 



D{t- M^ = pD{T\Mpre > M^; M,) + qpD{T\Mpre > M^; AQ 

* D{T\Mpre < M^; Me) + q^pD{T\Mpre > M'^; Me) 

* D{T\Mpre < M'^; Me) * D{T\Mpre < M^; Me) + ■ ■ ■ 

= pD{T\Mpre > M'^- Me) * ^Zo q^^lD^^Mpre < M^; Me)Y 



(14) 



where * denotes the convolution product and p is the probability that an 
earthquake is above M^, given that it is above Mc, i.e., 

p = Prob[A/ > M'^\M > Me] = iq-Hm',"^^-) ^ (15) 

using the GR law, whereas q = Prob[Af < M^\M > Me] = 1-p. Equation lO 
enumerates the number of ways in which two consecutive events for M > M^ 
may be separated by a recurrence time r; these are the number of events with 
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M < M'^ in between, each one contributing with a probability q, and then 
the time r between the two events is in fact a (fc + 1)— th return-time for the 
process with M > M^; from here and the independence between recurrence 
times the convolutions arise. 

Let us translate Eq. ffT^ to Laplace space, by using F{s) = Jg°° e^''^F{T)dT; 
then, the convolutions turn out to be simple products, i.e., 

oo 

D{s; AQ = pD{s\Mpre > Af^; A'Q ^ q''[D{s\Mpre < M'^; Me)]'=. (16) 

fe=0 

As qD{s\Mpre < M^; Mc) < 1 the series can be summed, yielding 
D(s- M') = pD{s\Mp.e>M-M.) 

using that D{s;Mc) = pD{s\Mpre > M^;Mc) + qD{s\Mpre < M'^]Mc). We 
have obtained an equation for the transformation of the recurrence-time prob- 
ability density under the thinning or decimation caused by the raising of the 
magnitude threshold from Mc to M'^. The second part in the process is the 
simple rescaling of the distributions, to make them have the same mean and 
comparable with each other; we obtain this by removing the effect of the 
decreasing of the rate, which, due to thinning, is proportional to p, so, 

D{T-M'^)^p-^D{p-\-NQ, (18) 

and in Laplace space, 

D{s-M'c)^D{ps-M'c). (19) 

Finally, the renormalization-group transformation T is obtained by combining 
the decimation with the scale transformation, 

T[Di^s- Mc)] = pD[ps\Mp.c^M'c-Mc) ^_ 

^ ^ '^ l-D{ps;Mc)+pD{ps\Mpre>M^;Mc) ^ ' 

A third step which is usual in renormalization-group transformations is the 
renormalization of the field, M in this case, but as we are only interested 
in recurrence times it will not be necessary here. The fixed points of the 
renormalization-group transformation are obtained by the solutions of the 
fixed-point equation 

T[D{s-Mc)]=D{s-Mc). (21) 

This equation is equivalent to the scaling law for the recurrence-time densities, 
Eq. 0), the only difference is that now it is expressed in Laplace space, as we 
are not able to provide the form of the operator T in real space. 
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4.2 Renormalization-Group Invariance of the Poisson Process 

We can get some understanding of the transformation T by considering first 
the simplest possible case, that in which there are no correlations in the pro- 
cess; so we have to break the statistical dependence between the magnitude 
and the subsequent recurrence time. This means that 

D{r\Mpre > M^; M,) = D{T\Mj,r.e < M^; M,) = D{t; M,) = Doir; A'Q 

(22) 
and then the renormalization transformation turns out to be 

-rrn f ^,t m pDa{ps]Mc) . 

^ms,M^)]-^_^j,^^^,.M^y (23) 

if we introduce lu = ps and substitute p — uj/s and q = 1 — lu/s in the 
fixed-point equation TDq{s;Mc) = Dq{s\Mc), we get, separating variables 
and equaling to an arbitrary constant k 

^ r - - = ^ r -- = k; (24) 

due to the fact that p and s are independent variables and so are s and lo. 
The solution is then 

Do{s;A'Q = {l + ks)-\ (25) 

which is the Laplace transform of an exponential distribution, 

Do{T;M,)^k-'e-^/\ (26) 

The dependence on Mc enters by means of k, a,s k = {t{Mc)); in the case of 
seismicity the GR law holds and k = R'^iMc) = Rgho'''^-. 

Summarizing, we have shown that the only process without correlations 
which is invariant under a renormalization-group transformation of the kind 
we are dealing with is the Poisson process. This means that if one consid- 
ers as a model of seismicity a renewal process (i.e., independent identically 
distributed return times) with uncorrelated magnitudes, then the recurrence- 
time distributions will not verify a scaling law when the threshold Mc is raised, 
except if D(t] Mc) is an exponential (which constitutes the trivial case of a 
Poisson process). 

Even further, the Poisson process is not only a fixed point of the transfor- 
mation, but a stable one (or attractor) for a thinning transformation in which 
events are randomly removed from the process (random thinning). If mag- 
nitudes are assigned to any event independently of any other variable (other 
magnitudes or recurrence times) the decimation of events after the risen of the 
threshold Mc is equivalent to a random thinning, and therefore the resulting 
process must converge to a Poisson process, under certain conditions |27j. 

The fact that for real seismicity the scaling function / is not an exponen- 
tial tells us that our renormalization-group transformation is not performing 
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a random thinning; this means that the magnitudes are not assigned inde- 
pendently on the rest of the process and therefore there exists correlations in 
seismicity. This of course is not new, but let us stress that correlations are 
fundamental for the existence of the scaling law Q : the only way to depart 
from the trivial Poisson process is to consider correlations between recur- 
rence times and magnitudes in the process. This is the motivation for the 
model explained in this section, for which we have chosen the simplest form 
of correlations between magnitudes and subsequent recurrence times. In fact, 
Molchan has shown that even for this correlated model the Poisson process is 
the only possible fixed point, implying that this type of correlations are too 
weak and one needs a stronger dependence of the recurrence times on history 
to depart from the Poisson case. After all, this is not surprising, as we know 
from the study of equilibrium critical phenomena that in order to flow away 
from trivial fixed points, long-range correlations are necessary. Therefore, the 
problem of finding a model of correlations in seismicity yielding a nontrivial 
recurrence-time scaling law is open. 



5 Correlations in Seismicity 

In the preceding section we have argued that the existence of a scaling law 
for recurrence time distributions is inextricably linked with the existence of 
correlations in the process, in such a way that correlations determine the form 
of the recurrence-time distribution. In consequence, an in-depth investigation 
of correlations in seismicity is necessary. 

Our analysis will be based in the conditional probability density; for in- 
stance, for the recurrence time we have, 

Prob[T < recurrence time < t + dr \ X] 
D{r\X) ^ , 

where \X means that the probability is only computed for the cases in which 
the condition X is fulfilled. If it turns out to be that D{t\X) is undistin- 
guishable from the unconditional density, D{t), then, the recurrence time is 
independent on the condition X; on the contrary, if both distributions turn 
out to be significantly different, this means that the recurrence time depends 
on the condition X and we could define a correlation coefficient to account 
for this dependence, although in general we might be dealing with a nonlinear 
correlation. 

Moreover, as we will compare values of the variables in different times 
(for example, the dependence of the recurrence time r^ on the value of the 
preceding recurrence time, Ti_i, for all i), we introduce a slight modification 
in the notation, particularly with respect the previous section, including the 
subindiccs denoting the ordering of the events in the probability distributions. 
Further, in order to avoid complications in the notation, we will drop the 
dependence of the conditional density on M^ when unnecessary. 
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5.1 Correlations between recurrence times 

Let us start with the temporal sequence of occurrences, for which we obtain 
the conditional distributions D{Ti\Ta < t^^i < Tf,); in particular we distin- 
guish two cases; short preceding recurrence times, D{Ti\Ti-i < Tb), where Tb 
is small, and long preceding recurrences, D(ri|Ti_i > Ta), with Ta large. The 
results, both for worldwide seismicity and for Southern-California stationary 
seismicity, turn out to be practically the same, see Fig. 6 of Ref. ^J. For 
short Ti-i, a relative increase in the number of short r^ and a decrease of long 
Tj is obtained, in comparison with the unconditional distribution, which leads 
to a steeper power-law decay of the conditional density for short and interme- 
diate times. In the opposite case, long ri_i's imply a decrease in the number 
of short Ti and an increase in the longer ones, in such a way that a flatter 
power-law exists here. In any case, the behavior for long r^ is exponential. So, 
short Ti-i's imply an average reduction of n and the opposite for long r^-i's, 
and then both variables are positively correlated. 

This behavior corresponds to a clustering of events, in which short re- 
currence times tend to be close to each other, forming clusters of events, 
while longer times tend also to be next each other. This clustering effect is 
different from the clustering reported in previous sections, associated to the 
non-exponential nature of the recurrence-time distribution, but is similar, in 
some sense, to the usual clustering of aftershock sequences, as these sequences 
also show this kind of correlations, although mainly due to the time-variable 
rate. 

In fact, the case of nonstationary seismicity was studied by Livina et al. 
for Southern California |S21 OHji with the same qualitative behavior. These 
authors explain their results in terms of the persistence of the recurrence time, 
which is a concept equivalent to the kind of clustering we have described. The 
results could be also similar to the long-term persistence observed in climate 
records. 

The effect of correlations can be described in terms of a scaling law, 
which constitutes a generalization of the scaling law for (unconditioned) re- 
currence time distributions, Eq. Q). In this way, we can write the conditional 
recurrence-time distribution in terms of a scaling function which depends on 
two variables, Rri and Rra, or Rri and Rrb^ see Ref. 33 . Further, the study 
of correlations between recurrence times can be studied beyond consecutive 
events, i.e., we can measure the distribution of Ti conditioned to Ti^2, or T^^a, 
etc. The results for these distributions show no qualitative difference, at least 
up to I — 10, in comparison with what we have explained for i — 1. 

The main results of this subsection can be summarized in one single sen- 
tence, reflecting the positive correlation between recurrence times: the shortest 
the time between the two last earthquakes, the shortest the recurrence of the 
next one, on average. 
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5.2 Correlations between recurrence time and magnitude 

If magnitudes are taken into account, there are two main types of correlations 
with the recurrence times. First, we consider how the magnitude of one event 
influences the recurrence time of a future event, in particular the next one, 
measuring D{Ti\Mi-i > M^). The results for the case of worldwide seismicity 
and for Southern California (in a stationary case) are again similar and show 
a clear (negative) correlation between Ali-i and Ti |H4| . 

Figure deleft) shows, for Southern-California in the stationary case, how 
larger values of the preceding magnitudes, given by Mi_i > M^, lead to a 
relative increase in the number of short r^ and a decrease in long Ti , implying 
that Mi-i and r^ are anticorrelated. For the cases for which the statistics 
is better, the densities show the behavior typical of the gamma distribution, 
with a power law that becomes steeper for larger M^; the different values of 
the power-law exponent are given at the figure caption. 
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Fig. 7. (left) Recurrence-time distributions for Southern-California conditioned to 
the value of the preceding magnitude, for the stationary period comprised between 
1994-1999. Each set of curves, which have been shifted up for clarity sake, corre- 
sponds to a value of Mc, which is, from bottom to top: 1.5, 2, 2.5, 3, 3.5. (right) 
Rescaled recurrence-time distributions conditioned to the value of the preceding 
magnitude, both for worldwide and for Southern-California stationary seismicity. 
Each set of data, shifted up again for clarity sake, corresponds this time to a differ- 
ent value of M^ — Mc, these being 0, 0.5, 1, 1.5, and 2, from bottom to top, and is 
fit by a gamma distribution, with (decreasing) power-law exponents 0.30, 0.45, 0.52, 
0.65, and 0.77, respectively. 



Remarkably, this behavior can be described by a scaling law for which 
the scaling function depends now on the difference between the threshold 
magnitude for the i — 1 event, M^, and the threshold for i, M^, i.e.. 



with R{Mc, M^) 



= l/(r(Mc,Af^)) and {Ti{Mc,M^)} the mean of the distri- 
bution D{Ti\AIi^i > M^; Mc). The original scaling law for unconditioned dis- 
tributions, Eq. Q, is recovered taking the case M^ = Mc- Figure [Tf right) 
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illustrates this scaling law both for worldwide and for Southern-California 
stationary seismicity. 

The second kind of correlations deals with how the recurrence time to 
one event r^ influences its magnitude Mi, or, equivalently, how the recurrence 
time to one event depends on the magnitude of this event. For this purpose, 
we measure D{Ti\Mi > M'^); the results for Southern-California stationary 
seismicity are shown in Fig.|HIleft). It is clear that in most of their range the 
distributions are nearly identical, and when some difference is present this is 
inside the uncertainty given by the error bars. However, there is one exception: 
very short times, r^ ~ 3 min, seem to be favored by larger magnitudes, or, 
in other words, short times lead to larger events; nevertheless, due to the 
short value of the time involve, we can ignore this effect. Then, the recurrence 
time and the magnitude after it can be considered as independent from a 
statistical point of view, with our present resolution (it might be that a very 
weak dependence is hidden in the error bars of the distributions) . 

Further, as in the previous subsection, we have gone several more steps 
backwards in time, measuring conditional distributions up to D{Ti\Mi^io > 
M^), and also we have extended the conditional distributions to the future, 
measuring D{Ti\Mj > M^) with j > i. The results are not qualitatively dif- 
ferent than what is described previously, with the first type of distributions 
dependent on M^ and the second type independent. From here we can con- 
clude the dependence of the recurrence time and the independence of the 
magnitude with the sequence of previous recurrence times, at least with our 
present statistics and resolution. 

Two sentences may serve to summarize the behavior of seismicity described 
here. First, the bigger the size of an earthquake, the shortest the time till next, 
due to the anticorrelation between magnitudes and forward recurrence times. 
Note that in the case of stationary seismicity this result is not trivially derived 
from the law of aftershock productivity. Second, the belief that the longer the 
recurrence time for an earthquake, the bigger its size, is false, as the magnitude 
is uncorrelated with the previous recurrence times. This shows clearly the time 
irreversibility of seismicity. 

5.3 Correlations between magnitudes 

Although not directly related with the temporal properties, we study the 
correlations between consecutive magnitudes, Afi_i and Mi, by means of the 
distribution D(Mi\Mi_i > M'^). As the analysis of the 1994-1999 period for 
Southern California did not provide enough statistics for the largest events, 
we considered a set of stationary periods, these being: Jan 1, 1984 - Oct 15, 
1984; Oct 15, 1986 - Oct 15, 1987; Jan 1, 1988 - Mar 15, 1992; Mar 15, 1994 
- Sep 15, 1999; and Jul 1, 2000 - Jul 1, 2001; all of them visible in Fig. El 

Again we find that both worldwide seismicity and stationary Southern- 
California seismicity share the same properties, but with a divergence for short 
times. Figure 10 of Ref. ^I] shows the distributions corresponding to the two 
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Fig. 8. (left) Recurrence-time distributions for Southern-California conditioned to 
the value of the magnitude of the incoming event, for the stationary period comprised 
between 1994-1999. The curves have been shifted up and each set corresponds to a 
value of Mc, which ranges from 1.5 to 3.5 in steps of 0.5, from bottom to top. (right) 
Distribution of magnitudes conditioned to the value of the preceding magnitude for 
Southern California for the set of stationary periods indicated at the text. Only 
events separated by r > 30 min have been considered. 



regions, and whereas for the worldwide case the differences in the distributions 
for different M^ are compatible with their error bars, for the California case 
there is a systematic deviation, implying a possible correlation. 

In order to find the origin of this discrepancy we include an extra con- 
dition, which is to restrict the events to the case of large enough recurrence 
times, so we impose Ti > 30 min. In this case, the differences in Californian 
distributions become no significant, see Fig. ISJright), which means that the 
significant correlations between consecutive magnitudes are restricted to short 
recurrence times [SI]. Therefore, we conclude that the Gutenberg-Richter law 
is valid independently of the value of the preceding magnitude, provided that 
short times are not considered. This is in agreement with the usual assumption 
in the ETAS model, in which magnitudes are generated from the Gutenberg- 
Richter distribution with total independence of the rest of the process |18II19| . 
Of course, this independence is estabhshed within the errors associated to our 
finite sample. It could be that the dependence between the magnitudes is 
weak enough for that the changes in distribution are not larger than the un- 
certainty. With our analysis, only a much larger data set could unmask this 
hypothetical dependence. 

The deviations for short times may be an artifact due to the incomplete- 
ness of earthquake catalogs at short time scales, for which small events are 
not recorded. Helmstetter et al. |3S] propose a formula for the magnitude of 
completeness in Southern California as a function of the elapsed time since a 
mainshock and its magnitude; applying it to our stationary periods, for which 
the larger earthquakes have magnitudes ranging from 5 to 6, we obtain that 
a time of about 5 hours is necessary in order that the magnitude of com- 
pleteness reaches a value below 2 after a mainshock of magnitude 6. After 
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mainshocks of magnitude 5.5 and 5 this time reduces to about 1 hour and 
15 min, extrapolating Helmstetter et al.'s results. In any case, it is perfectly 
possible that large mainshocks (not necessarily the preceding event) induce 
the loss of small events in the record and are the responsible of the deviations 
from the Gutenberg-Richter law at small magnitudes for short times. If an 
additional physical mechanism is behind this behavior, this is a question that 
cannot be answered with this kind of analysis. 

As in the previous subsection, we have performed measurements of the 
conditional distributions for worldwide seismicity involving different Mi and 
Mj, separated up to 10 events, with no significant variations in the distribu- 
tions, as expected. This confirms the independence of the magnitude Mi with 
its own history. 

These results, together with those of the previous subsection allow to state 
that an earthquake does not "know" how big is going to he (at least from the 
information recorded at the catalogs, disregarding spatial structure, and with 
our present resolution) [5HITT] . 

5.4 Correlations between recurrence times and distances 

Up to now we have considered seismicity as a point process in time, marked by 
the magnitude. If, in addition to this, we take into account the spatial degrees 
of freedom, these new variables allow to study other types of correlations. 
Of outstanding importance will be the distances between earthquakes, and 
specially the distances between consecutive earthquakes, which we may call 
jumps. 

The correlations between jumps and recurrence times have been investi- 
gated in Ref. {36! !, and they show a curious behavior. There are two kinds of 
recurrence time distributions conditioned to the distance, one for short dis- 
tances, which can be represented by a gamma distribution with a decaying 
power law of exponent around 0.8, and the distribution for long distances, 
which is an exponential. It is clear then that in one case we are dealing with 
aftershocks and in the other with Poissonian events. The particularity of these 
distributions is that they are independent on the distances, provided that the 
set of values of the distances are short, or long. For worldwide earthquakes, the 
difference between short and long distances is around 2° (200 km), whereas 
for Southern California this value is 0.1°, approximately. 

We may note that the (unconditional) distribution of recurrence times is 
then a mixture of these two kind of conditional distributions, and therefore, 
the existence of a universal recurrence-time distribution is a consequence of a 
constant proportion of short and long distances in seismicity, or of aftershocks 
and uncorrelated events. 
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6 Bak et al.'s Unified Scafing Law 

Bak, Christensen, Danon, and Scanlon introduced a different way to study 
recurrence times in eartliquakes ^. Tliey divided tlie region of Southern Cal- 
ifornia into approximately equally-sized squared subregions (when longitude 
and latitude are taken as rectangular coordinates), and computed the series of 
recurrence times for each subregion. The main difference with the procedure 
explained in the previous sections is that Bak et al. included all the series of 
recurrence times into a unique recurrence-time distribution, performing there- 
fore a mixing of the distributions for all subregions. As seismic rate displays 
large variations in space (compare the rates of occurrence in Tokyo and in 
Moscow, and the same happens at smaller scales) Bak et al.'s procedure leads 
to a very broad distribution of recurrence times. 

It was found that the recurrence-time densities defined in this way, 
'D{t;Mc,£), for different magnitude thresholds Mc and different linear size 
i of the subregions, verify the following scaling law (see Fig. I^Ileft)), 

V{T;M,,e)^TZF{nT), 

which was named unified scaling law, where F is the scaling function, showing 
a power-law decay with exponent close to 1 for small recurrence times and 
a different power-law decay for long times, with exponent around 2.2 [37], 
whereas TZ{Mc,£) is the spatial average of the mean seismic rate, i.e., the 
average of Rxy{Mc, t) for all the regions with seismic activity (labeled by xy), 
so, TZ = ^^ Rxy/n, where n is the number of such regions. Note that TZ 
is the inverse of the mean of V. From the GR law for each region, R^y — 
Rxyo^^'^^'" and from the fractal scaling of n with i, n = {C/iYf , we get, 
TZ = Roil/CynO-''^'''', with Ro{C) = Y^xyR^yoi^) and C a rough measure 
of the linear size of the total area under study. Therefore we can write the 
scaling law as 

which relates the recurrence-time density, defined in the Bak et al.'s way, with 
the GR law and with the fractal distribution of epicenters, and from here the 
name of unified scaling law. Molchan and Kronrod have studied this law in 
the framework of multifractals 38 . 

Later it was found that the unified scaling law holds beyond the case of 
Southern California, for instance for Japan, Spain, New Zealand, New Madrid 
(USA), or Iceland, as well as worldwide 39, 40 . However, it turned out that 
the scaling function is not universal, as there are differences for different re- 
gions, mainly in the crossover between short and long times, although the 
value of the long-time power-law exponent seems to be in all cases 2.2, and 
therefore universal, see Fig. Efleft) . The deviations from the hyperbolic- like 
behavior (exponent close to one) for very short times have also been studied 
[40. 

It is clear that the short-time exponent must be related (but not identi- 
cal!) to the Omori p— value; on the other hand, the long-time exponent is a 
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Fig. 9. (left) Mixed recurrence-time densities (defined in the way of Bak et al.) and 
rescaled by TZ, for several values oil and Mc- The different sets of curves correspond, 
from top to bottom, to (i) Southern California, 1984-2001; (ii) Northern California, 
1985-2003; (iii) Stationary seismicity: Southern California, 1988-1991; worldwide, 
1973-2002; Japan, 1995-1998; and Spain, 1993-1997; (iv) Stationary seismicity: New 
Zealand, 1996-2001, and New Madrid, 1975-2002; A total of 84 distributions are 
shown, I ranging from 0.039° to 45°, and 1.5 < Mc < 6. The distributions are shifted 
to the bottom for clarity sake. All the left tails are fit by a (decreasing) power-law 
with exponent 2.2, the right part of the distributions are fit by a power law with 
exponent 0.95 or 0.9, see Ref. |39|. (right) Distribution of mean rates of occurrence 
for events with M > 2 in Southern-CaUfornia, 1984-2001 (averaging 1984-1992 and 
1993-2001), using diverse values of (.. The distributions are rescaled by L^'^ and the 
straight line is a power law fit which turns out to be oc 1/ BP'^ . 



consequence of a power law distribution of seismic rates in space, as we now 
show 131] ■ Therefore, with the purpose of understanding the relation of Bak et 
al.'s results with the rest of this work, let us generalize the nonhomogeneous 
Poisson-Omori process previously introduced, in order to include the univer- 
sal scaling law for Omori sequences. We have explained that these sequences 
can be characterized by an r— dependent recurrence-time probability density 
of the form D{T\r) oc r'''r'''~^e^'"^/° (note that this includes the nonhomoge- 
neous Poisson process, given by 7 = a = 1, but for real field data 7 ~ 0.7). We 
expect that, for a given spatial area, this is valid not only for Omori sequences 
but also for a general time- varying rate; then the overall probability density of 
the recurrence times, independently of r, is given by the mixing of all D{T\r) 



1 /"■" 

D{T\r^) = - / rD 



{T\r)p{r)dr, 



(27) 



in fact, this is just Eq. ©; we recall that p{r) is the density of rates, p is the 
mean rate, rg is the maximum rate, and r^ is the minimum rate, related to the 
background seismicity level. Note that we have emphasized the dependence 
of the resulting distribution on r„j. 

Let us consider that the distribution of rates comes essentially from Omori 
sequences, then, as we already know, p{r) = C/r^+^Z^. The analysis is sim- 
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plified for 7 = 1/p, although the conclusions will be of general validity, so, in 
this case, 

D{r\r^) oc - ^ -.^^ 1, (28) 

where, in the same way as for a nonhomegeneous Poisson process, the mini- 
mum rate r^ determines the exponential tail of D{T\rm) for large r, which is 
preceded by a decreasing power law with exponent 2 — 1/p ii ro ^ rm- 

Up to now we have arrived to a slightly different, more convenient varia- 
tion of the distribution corresponding to a nonhomogeneous Poisson process. 
Next step is to take into account the spatial degrees of freedom, fundamen- 
tal in Bak et al.'s approach. In fact, as we have explained, their approach 
performs a mixing of recurrence times coming from different spatial areas (or 
subregions), which are characterized by disparate seismic rates. In particular, 
each area will have a different r^, depending on its background seismicity 
level. As the minimum rate is difficult to measure (it depends on the size of 
the time intervals selected) , we assume that the minimum rate r™ is somehow 
proportional to the mean rate of the sequence fi, which in turn is in corre- 
spondence with the mean rate in the area, R. This spatial heterogeneity of 
seismicity can be well described by a power-law probability density of mean 
rates R, p{R) oc l/i^^-", with a ~ 0.2, see Fig. EIright) and Ref. [37|; then, 

p(r,„) oc l/r,V" 
and therefore the recurrence-time probability density comes from the mixing. 



2?(t) oc / rmD{T\rm)p{rm)drm (29) 

where r^ varies between rmm and rmM- Integration, taking into account that 
C//J, depends on r„, leads, for rmmT ^ 1 ^ t^mt, to 

V(t) oc l/r2+" (30) 

In this way the power law for long times, reflects the spatial distribution 
of rates. The universal value of the exponent 2 + a ^^, would imply the 
universality of seismicity spatial heterogeneities. In consequence, Bak et al.'s 
unified scaling law provides a way to measure these properties. Further, Eq. 
H30|l shows that the change of exponent in 2?(t) appears for r larger than 
l/r„iM, which corresponds, for the area of highest seismicity, to the mean of 
events that are in the tail of the Omori sequence, or in background seismicity, 
and therefore at the onset of correlation with the mainshock. It is in this sense 
that the change of exponent separates events with different correlation. On 
the other hand, the power law for short times is not affected by the spatial 
mixing and therefore V oc l/r^"^/^. 
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7 Conclusions 

We hope we have convinced the reader about the interest to study of the 
temporal features of seismicity. The research of the author was illuminated by 
the ideas and philosophy of the late Per Bak. 
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